% By Martin Andreasen, 2009
% This function computes second order analytical derivatives 
% of the pricing kernel. The integer "R_function" has the following properties:
% 1  - if log-transformation of bond prices is used and 
% 0  - if no log-transformation of bond prices is used

function [M,Mx,Mxp,My,Myp,...
    Mxx,Mxxp,Mxy,Mxyp,Mxpx,Mxpxp,Mxpy,Mxpyp,Myx,Myxp,Myy,Myyp,Mypx,Mypxp,Mypy,Mypyp] = ...
    Anal_PricingKernel_Real_derivatives(x,xp,y,yp,order_app,R_function)

% Some integers
nx    = size(x,2);
ny    = size(y,2);

%Define parameters
syms BETTA B CHI CHI0 THETA DELTA ALFA PHI PHIzero ZI KAPAw ETA RHOR PHIpai PHIpai_1 ...
     PHIy PHIy_1 PHIc PHIc_1 PHIl NU U0 U0d...
     RHOz  RHOd   RHOn  RHOp RHOa ...
     Kss OUTPUTss PAIss MUZss Css AA lss Wss;
Rss = sym('Rss');
symparams=[BETTA,B,CHI,CHI0,THETA,DELTA,ALFA,PHI,PHIzero,ZI,KAPAw,ETA,...
     RHOR,PHIpai,PHIpai_1,PHIy,PHIy_1,PHIc,PHIc_1,PHIl,NU,U0,U0d,...
     RHOz ,RHOd   RHOn  ,RHOp,RHOa, ...
     Kss OUTPUTss PAIss MUZss Css AA lss,Rss,Wss];

%Define variables 
syms c_cu   r_cu    pai_cu   w_cu  l_cu  output_cu  mc_cu  evf_cu  vf_cu  p1_cu  p1Real_cu  dc_cu...
     c_cup  r_cup   pai_cup  w_cup l_cup output_cup mc_cup evf_cup vf_cup p1_cup p1Real_cup dc_cup...
     c_ba1  muz_cu   d_cu  n_cu  paistar_cu  a_cu...
     c_ba1p muz_cup  d_cup n_cup paistar_cup a_cup;
  
% The real stochastic discount factor
mReal_cup = BETTA*d_cup/d_cu*...
    (AA*evf_cu^(1/(1-ALFA))/(vf_cup*muz_cup^((1-CHI)*(1-CHI0)) ))^ALFA...   
    *(c_cup-B*c_cu/muz_cup)^(-CHI)/((c_cu-B*c_ba1/muz_cu)^(-CHI))*muz_cup^(-CHI*(1-CHI0)-CHI0);  


% The pricing kernel
M = mReal_cup;

%Make f a function of the logarithm of the state and control vector
if R_function == 1
    y_log  = [c_cu   r_cu    pai_cu   w_cu  l_cu  output_cu  mc_cu  evf_cu  vf_cu];
    yp_log = [c_cup  r_cup   pai_cup  w_cup l_cup output_cup mc_cup evf_cup vf_cup];
    x_log  = [c_ba1  muz_cu   d_cu  n_cu  paistar_cu  a_cu];
    xp_log = [c_ba1p muz_cup  d_cup n_cup paistar_cup a_cup];

    M = subs(M, [x_log,y_log,xp_log,yp_log], exp([x_log,y_log,xp_log,yp_log]));
end
% For Extended Path
M_ep = char(M);
M_ep = strrep(M_ep,'*','.*');
M_ep = strrep(M_ep,'/','./');
M_ep = strrep(M_ep,'^','.^');

% The first order derivatives of the pricing kernel
% Recall we only need the first order derivatives of the pricing kernel to get
% bond prices up to second order
if order_app > 1
    Mx  = jacobian(M,x);
    Mxp = jacobian(M,xp);
    My  = jacobian(M,y);
    Myp = jacobian(M,yp);
else
    Mx  = zeros(1,nx);
    Mxp = zeros(1,nx);
    My  = zeros(1,ny);
    Myp = zeros(1,ny);
end

% The second order derivatives of the pricing kernel
% Recall we only need the second order derivatives of the pricing kernel to get
% bond prices up to third order
if order_app > 2 
    % With respect to Mx
    Mxx   = reshape(jacobian(Mx(:) ,x) ,nx,nx);
    Mxxp  = reshape(jacobian(Mx(:) ,xp),nx,nx);
    Mxy   = reshape(jacobian(Mx(:) ,y) ,nx,ny);
    Mxyp  = reshape(jacobian(Mx(:) ,yp),nx,ny);

    % With respect to Mxp
    Mxpx  = reshape(jacobian(Mxp(:),x) ,nx,nx);
    Mxpxp = reshape(jacobian(Mxp(:),xp),nx,nx);
    Mxpy  = reshape(jacobian(Mxp(:),y) ,nx,ny);
    Mxpyp = reshape(jacobian(Mxp(:),yp),nx,ny);

    % With respect to My
    Myx   = reshape(jacobian(My(:) ,x) ,ny,nx);
    Myxp  = reshape(jacobian(My(:) ,xp),ny,nx);
    Myy   = reshape(jacobian(My(:) ,y) ,ny,ny);
    Myyp  = reshape(jacobian(My(:) ,yp),ny,ny);

    % With respect to Myp
    Mypx  = reshape(jacobian(Myp(:),x ),ny,nx);    
    Mypxp = reshape(jacobian(Myp(:),xp),ny,nx);
    Mypy  = reshape(jacobian(Myp(:),y),ny,ny);
    Mypyp = reshape(jacobian(Myp(:),yp),ny,ny);
else
    % With respect to Mx
    Mxx   = zeros(nx,nx);
    Mxxp  = zeros(nx,nx);
    Mxy   = zeros(nx,ny);
    Mxyp  = zeros(nx,ny);

    % With respect to Mxp
    Mxpx  = zeros(nx,nx);
    Mxpxp = zeros(nx,nx);
    Mxpy  = zeros(nx,ny);
    Mxpyp = zeros(nx,ny);

    % With respect to My
    Myx   = zeros(ny,nx);
    Myxp  = zeros(ny,nx);
    Myy   = zeros(ny,ny);
    Myyp  = zeros(ny,ny);

    % With respect to Myp
    Mypx  = zeros(ny,nx);    
    Mypxp = zeros(ny,nx);    
    Mypy  = zeros(ny,ny);    
    Mypyp = zeros(ny,ny);
end
